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Abstract 

In the analysis of vibrations of continuous elastic systems, one often encounters 
complicated transcendental equations with roots directly related to the system’s natural 
frequencies. Typically, these equations contain system parameters whose values must be 
specified before a numerical solution can be obtained. The present paper presents a 
method whereby the fundamental frequency can be obtained in analytical form to any 
desired degree of accuracy. The method is based upon truncation of rapidly converging 
series involving inverse powers of the system natural frequencies. A straightforward 
method to developing these series and summing them in closed form is presented. It is 
demonstrated how Computer Algebra can be exploited to perform the intricate analytical 
procedures which otherwise would render the technique difficult to apply in practice. We 
illustrate the method by developing two analytical approximations to the fundamental 
frequency of a vibrating cantilever carrying a rigid tip body. The results are compared to 
the numerical solution of the exact (transcendental) frequency equation over a range of 
system parameters. 
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Introduction 

The general availability of computer algebra systems has resulted in analyses of 
complicated problems which heretofore have been regarded as analytically intractable.These 
tools practically eliminate the tedious error prone manipulations required by hand- 
derivations and allow the analyst to explore various analytical treatments which would be 
too costly otherwise. In the same way that digital simulation has revolutionized the 
numerical treatment of engineering problems, symbolic computation promises to be a 
powerful tool in analytical investigations. 

In the area of multibody dynamics, computer algebra has been used to derive the full 
nonlinear equations of motion in symbolic form [1]. These equations are typically so 
complex that they daunt inspection. However, built in translators convert these equations 
into a higher programming language e.g FORTRAN which results in extremely efficient 
digital simulations. For sufficiently simple systems, symbolic representations can be of 
direct use in studying such issues as elastic stability and buckling [2], Perturbation methods 
are ideally suited to treatment by symbolic computation [3]. 

The classical method of determining the natural frequencies of a continuous elastic 
system results in an eigenvalue problem and associated characteristic equation. In general, 
this equation is transcendental and is embedded with various system parameters. Thus 
recourse must be made to numerical methods to solve these equations and the dependencies 
of the frequencies upon the various parameters can only be revealed through exhaustive 
computation. It therefore appears desirable to be able to approximate the roots of these 
equations by analytical expressions which are relatively simple yet accurate. The principal 
idea behind the method presented in this paper is to find closed form expressions for 
infinite series, the terms of which involve inverse powers of the natural frequencies. 
Truncating the series after the first term gives an approximation to the fundamental 
frequency. By summing sufficiently high powers, this approximation can be made 
arbitrarily accurate; but the resulting formula increases in complexity. The methods 
whereby others have addressed this problem are quite varied, ranging all the way from 
Fourier Series to complicated contour integration and difficult procedures involving integral 
equations. Hughes [4] obtains numerous modal identities by expanding the Green's 
function in a series of eigenfunctions. The technique we present is extremely simple and 
appears to have been applied in a restricted form by Lord Rayleigh [5]. Our result is very 
general and can be applied to any vibrational system one the characteristic equation is 
established. The only difficulty in applying the method is the need to develop complicated 
transcendental functions into Taylor series and manipulating the resulting coefficients. 
However, with the use of a computer algebra system, this task becomes almost trivial. The 
method is illustrated by deriving two analytical approximations to the fundamental 
frequency of a vibrating cantilever carrying a rigid tip body. The accuracy of these results is 
verified by comparisons with numerical solution of the frequency equation over a range of 
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parameter values. 


A pproximations Based Upon Rayleigh's Principle 

Before presenting the technique predicated on infinite series, let us consider a 
"symbolic" solution to a prototype vibration problem employing the celebrated Rayleigh 
Principle. The elastic system consists of an Euler-Bemoulli beam cantilevered at one end, 
and carrying a rigid tip body at the other (see Fig. l.A). The beam has a constant mass 
density (per unit length) p, uniform bending stiffness El , and nominal length Z. A rigid tip 
body of mass m and inertia J (about P) is attached to the beam tip at x=£. We denote by c 
the distance from P to the tip body mass center.The derivation of the partial differential 
equation of motion and associated eigenvalue problem is given in Appendix A. The system 
eigenvalues are the solutions Pk to eq.(A.9) and are seen to depend upon the three 
dimensionless tip body parameters 


rrT=J2-, r=-J- 

p£ 2 p£ 5 

The relationship between the eigenvalues and the system natural frequencies is given by 
eq.(A.lO). 

Let us approximate the beam deflection u(x,t) with a cubic polynomial in x. 


u(x,t)=£ 2 qi(t)+£ 3 q 2 (t) 


where t=\lZ and the geometric boundary conditions at x=0 have been observed. Here 
qi(t),q 2 (t) are undetermined generalized coordinates. The system kinetic energy T and 
potential energy V are then discretized into the respective quadratic forms (see eqs.(A.12) 
&(A.ll)) 


T = 


ei + m(l + 2c*) + i(l-m=2) 


qi + 


1 +3c +-SUj - me 2 ) 

14 2 2Z 2 


qi- 


m( 1 +2c*X 1 +3c*)+-^J-mc 2 ) 

6 Z 2 


qiQ2 


V = 2£/.( q 2 + 3 q ^ + 3qiq 2 ) 

Z 2 

If we write T=l/2q T [M]q and V= 1 /2q T [K] q , then the system’s first two natural 
frequencies are approximated by the roots (0; of the characteristic polynomial 
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det([K]-a) 2 [M])=0 


Expanding this determinant and solving the resultant quadratic is relatively painless if a 
computer algebra system is invoked.The resulting expression for the fundamental 


frequency can be written in the form 

oW^V 4 Pi 


with 


(3i = 


1260- 


1/4 


I(630c*+210)m*+630J*+t 3 +51. 
where t 3 =2Vl(33075 J* 2 + 1 260 J* +t 2 +ti +208) 1 h 
and t 2 =[(66150c*+l 1025)j*+4200c*+1680]m* 
t , =(44 1 00c* 2 + 22050c* + 3675)m* 2 


( 1 ) 


This result of course provides an upper bound to the true fundamental frequency. 

It should be noted that this method meets with practical difficulties when one attempts 
to improve the accuracy by retaining additional terms in the expansion of the elastic 
displacement. The higher degree of the concomitant characteristic polynomial renders an 
analytical solution impossible. The method to be described in the next section does not have 
this limitation. 


A pproximations Based Upon Infinite Series 

The current method is based upon truncation of infinite series in the frequencies (On 
such as^iT-i l W, l/w„ 4 etc. where the sum can be expressed as a relatively simple 
algebraic function of the system parameters. If the series convergence is sufficiently rapid, 
then truncating the series after the first term yields a formula which approximates the 
fundamental frequency coj. Clearly, by summing sufficiently high powers of (O^ 1 we can 
approximate the first frequency to any desired degree of accuracy and will always have a 
lower bound. As will be seen, the corresponding formulas become increasingly complex. 
However, it should be pointed out that the generation of these higher order results can 
always be carried out in practice unlike the procedure of the last section. Hughes [4] 
generates series like the above and expresses the sum as a volume integral containing 
products of the Green's function with the mass density. He notes the difficulty of 
performing these integrations when the powers of co^ 1 increase. Our method only requires a 
Taylor series expansion once the transcendental frequency equation is established. 

Before considering the case of a transcendental equation, we present an elementary 
result from the theory of polynomial equations. 


Given the polynomial equation 

l+otiz+a 2 z 2 +- • +a n z n =0 
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with roots z;(i=l ,2,- • •, n) (over the field of complex numbers), we can show that 


(a) £i=- a i 

i=l Zl 
n 

(b) S^=a?- 2a 2 
i=i zf 


Proof 

First note that if 0 is a root, it can be removed leaving a deflated polynomial with no 
zero roots. Hence there is no loss in generality if we assume all z;*0. The general 
polynomial with roots Zi,z 2 ,---,Zn can be written as 

(z-zi)(z-z 2 )- • (z-zn)=0 


Expanding and dividing through by the product (- l) n ziz 2 - • z,, we obtain 


1-(_L+_L+- . .+X| Z+ /_1_+_L_ + , . .+ — 1 — ] z 2+. . . + — — z n - o 

\Zl Z 2 Zfi ! \ZjZ 2 ZjZ3 Zj|_|Zf| / Zl Z 2 * ■ * Zn 


which is of the desired form. It follows that 


-aj =sum of reciprocals of roots 

a 2 =sum of products of the reciprocals of the roots taken 2 at a time 

n 

Hence, ai~2a 2 =^ . 

i=i A 

Sums involving higher powers of the inverse roots can also be generated. Thus 

X -V=3aia 2 -3a 3 -a5, 

i=i A 

n 

^ -L=ai-4aia 2 +2a2+4aia3-4a4 etc. 

i=i A 

In an effort to adopt these results to the case of a transcendental equation f(z)=0 with an 
infinite number of roots, it is natural to expand f(z) in a power series and formally apply 
the above formulas to this "infinite degree" polynomial. It turns out that this procedure can 
be proven mathematically valid when f is an entire function of the complex variable z. We 
will not present the proof here but refer the reader to [6] for the necessary theory of entire 
(integral) functions. 


82 



As a means of illustration, we apply this technique to approximate the fundamental 
frequency of the beam-tipbody system considered above. Using power series expansions 
for the trigonometric and hyperbolic functions, the frequency equation (A.9) assumes the 
form 


1 --Ul2J‘+ 12mV+4rrT+l)p 4 + 


3^° m*+ 168)j*-420m* 2 c* 2 +56m*c*+8m*+ 



(2) 


Since the coefficients of p and p 2 are zero, we conclude that £ ~0 and £ -L=0. These 

n-1 n-1 

two results become immediately obvious, since, if Pk>0 is a root of eq.(A,9) then so are: 
-Pk> ‘Pk»-«Pk- In order to obtain series converging to a nonzero result, write eq.(2) as 

1 +aiP 4 +ct2P 8 +- • =0 


and form the auxiliary "polynomial" 


1 +<XiZ + (X2Z 2 +- • =0 (3) 

If Pk is a root of eq.(2), then Zk=Pk is a root of eq.(3). This artifice coalesces the quadruple 
of roots (Pk -Pk.»Pk ,-iPk} of eq.(2) into a single root of eq.(3). Applying our method to 
the auxiliary equation(3), we obtain 


and 



(4) 


X -L=((5880c* 2 +3360c*+560)m* 2 +[(l0080c*+2520)j*+728c*+264]m* + 

A (5) 

5040J* 2 + 504J* + 33 }/5040 


Truncating the above series after the first term leads to the respective approximations 
and 


12 \ 1/4 

12J'+12m'c'+4m'+l * 

(6) 

B,=( 5Q4Q l l/8 

lsi + s 2 +s 3 / 

(7) 


83 



where 

si =(5880c* 2 +3360c*+560)m*2 

s 2 =[( I0080c*+2520)J* +728c*+264]m* 

s 3 =504 J*( 1+ 10J*)+ 33 

Numerical Results 

Verification of the modal identities (4) & (5) is provided in Table 1 below. The 
eigenvalues Pn were generated by numerically solving the transcendental equation (A.9); 
the sequences of partial sums appear in the last two columns. The numbers in the last row 
(n=<») were obtained from the theoretical values appearing on the right hand sides of 
equations (4) and (5). All values were generated with m*=2.0, J*=0.028, and c*=0.1 . 


Table 1 - Partial Sums of Series: 

Eigenvalues of Beam/tipbody 


~n 

Pn 

■*951 

EH 

i 

1.0077 

.9699 

.9408 

2 

3.4599 

.9769 


3 

5.9100 

.9777 


4 

8.5047 

.9779 


5 

11.3806 

.9780 


oo 


.9780 
(eq. 4) 

BHKBlwl 


The two approximations to the "dimensionless frequency" Pi based upon series 
truncation (eqs. (6)&(7) ) are tabulated in Table 2 below for the case of a pure tip mass - 

J*=c*=0. The values in the second column (’Pf) were obtained by numerical solution of 
eq.(A.9). As the value of the tip mass increases (relative to the mass of the beam),both 
approximations improve. As expected, the approximation based upon eq.(7) is superior to 
that supplied by eq.(6). 
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Table 2 - Fundamental Frequency Approximations for Beam with Tip Mass 


m* 

Pi 

Eq. 6 (% error) 

Eq. 7 (% error) 

0.0 

3.5160 

3.4641 (1.5) 

3.5154 (1.6 E-02) 

5.0 

0.7569 

0.7559 (0.13) 

0.7569 (1.4 E-04) 

10.0 

0.5414 

0.5410 (.069) 

0.5414(3.7 E-05) 

15.0 

0.4437 

0.4435 (.047) 

0.4437 (1.7 E-05) 

20.0 

0.3850 

0.3849 (.035) 

0.3850 (9.7E-06) 


Table 3 below is similar in format to Table 2 but was generated with J*=l and c*=0, 
which represents a relatively large concentrated inertia at the tip of the beam. In this case we 
see that the approximations degrade with increasing m*. 

Table 3 - Fundamental Frequency Approximations for Beam with Tip Body 

c*=0, J*= 1 


m* 

Pi 

Eq. 6 (% error) 

Eq. 7 (% error) 

1.0 

.8679 

.8402 (3.2) 

.8669(.ll) 

1.3 

.8406 

.8120(3.4) 

.8396 (.12) 

1.5 

.8236 

.7947(3.5) 

.8225 (.13) 

1.8 

.7995 

.7708 (3.6) 

.7984 (.14) 

2.0 

.7845 

.7559 (3.6) 

.7833 (.14) 
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Appendix A 


We shall presently analyze the free vibration of a uniform cantilevered beam carrying a 
rigid tip body. Expressions for the potential and kinetic energies as well as the 
transcendental frequency equation are established. 

Fig. l.A below depicts the system in a deformed state. A uniform beam of mass 
density p (per unit length), bending stiffness El and length £ lies along the x axis when in 
equilibrium. A rigid body of mass m and moment of inertia J (about P) is attached to the 
beam tip at P. The distance between P and the rigid body mass center is c, and this 
directed line segment lies along the beam tip tangent direction (to prevent the tip body from 
exerting axial loads onto the beam). The small transverse displacement of the beam is 
denoted by u(x,t). 



We shall assume that the displacement of the beam and its slope are small quantities and 
therefore make the approximation that the angle 9p between the x axis and the beam tip 
tangent line at P can be approximated by u x (£,t). Denoting the inertial velocities of P and 
the tip body mass center by vp and v$ respectively, we can write 


V® = Vp+ CDpXC 


where cop=u xt (£,t)k is the angular velocity of the tip body and c is the vector from P to the 
tip body mass center. Recalling that |8p] is small and neglecting the term 9p0p, we find 


v© = 


^, 0+^,0 
9t dxdt 


(A.l) 


The expressions for the absolute translational and rotational accelerations of the tip body 


C-2 



mass center follow from the above by direct differentiation. 

In order to write the boundary conditions for u(x,t) at the endpoint x=£, we consider a 
free body diagram of the tip body. As indicated in Fig. 2.A, the beam exerts a force S 
directed along the y axis and a moment M directed along the z axis upon the tip body at 
the point P. 



Fig. 2.A free body diagram of tip body 


The equation of motion for the tip body along the y axis is 


S=m 




at 2 


3x3t 2 


From elementary beam theory the shearing force in the beam at x=£ is given by 

S=EI 9 3 u/ax 3 | x= £ . In conjunction with the above, this supplies one of the required 
boundary conditions. 


a 3 u 

El m| 

3x3 


3 2 u 


-+c- 


a 3 u 


3t 2 3x3t 2 


=0 at x=£ 


(A. 2) 


The second boundary condition at x=£ is obtained by considering the rotational motion of 
the tip body. If we denote by h the angular momentum of the tip body about its mass 
center, then we have the relation 

fh = M-cxS 
dt 

Taking the z component of this equation.neglecting the second order term in cop and using 
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the relation M=-EId' 2 u/9x 2 J x =£k, we arrive at the result 


-,2 -\2 -\3 

^dU du , du _ « 

El +mc — -+J =0 at x=*L 


dx 2 


dt 2 dxdt 2 


(A. 3) 


Since the beam is clamped at x=0, we have the two additional boundary conditions 

u(0,t)=0 and ^<0,t)=0 (A.4) 

dx 

The partial differential equation for free vibration is the well known relation 

3 4 -> 2 

_,d u du _ 

El + p =0 

9x 4 at 2 

We now proceed to solve the above homogeneous equation subject to the geometric 
boundary conditions (A.4) and natural boundary conditions (A.2) & (A.3). Seeking 
solutions of the form e^tpfx) we are led to the eigenvalue problem 


d4<p -^=0 

dx 4 

(A. 5) 

<p”(£) + ^|(p(£)+c(p(2)]=0 

(A. 6) 

<P*(£) - j^(mC <p(£) + J 9(2)] = 0 

(A.7) 

(p(0)=<p(0)=0 

(A.8) 

where (') indicates differentiation with respect to x. 


It can be shown that all the eigenvalues are positive. The general solution of eq.(A.5) is 


<p(x)=Cisin otx+C2Cos ax+C3sinhax+C4Coshax 

(a=X l/ 4 >o) 

In order to have a nontrivial solution satisfying the boundary conditions (A. 6), (A. 7) & 
(A.8), the eigenvalues must satisfy the transcendental characteristic equation 
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m*(j*-m*c* 2 )p 4 (l-cos Pcoshp) + m*p(cos P sinhP-sin P coshp) 
-2m*c*p 2 sin P sinhP-J*p 3 (sin P cosh P+sinhP cos p) (A.9) 

+ 1 +cos P coshp=0 

where we have introduced the "dimensionless frequency " P=a£ and the dimensionless tip 
body parameters are defined by 

m*= n y p £. c *=c/£, J*=J/p £ 3 

The natural frequencies are then given by 


COk = 


-EL pj 


(A. 10) 


For purposes of reference, the systems potential energy V is in the form of strain energy 
stored in the beam and is given by the formula 



Ul \2 
a u 

(ax 2 


dx 


(A. 11) 


while the kinetic energy T is the sum of the translational kinetic energy of the beam and tip 
body with the rotational kinetic energy of the tip body. Employing eq.(A.l) we can write 



(-[pdx+m 


Ul 2 

dt dxdt 


+±{j-mc 2 )| 


3 2 

o u 


dxdt 


< 2,0 


(A. 12) 
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